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£h ; Abstract 

i , 

^ A circular twist disclination is a nontrivial example of a defect in an elastic continuum that causes large de- 
q formations. The minimal potential energy and the corresponding displacement field is calculated by solving the 

' Euler-Lagrange-equations. The nonlinear incompressibility constraint is rigorously taken into account. 
I— ~~ ', By using an appropriate curvilinear coordinate system a finer resolution in the regions of large deformations is 
I ' obtained and the dimension of the arising nonlinear PDE's is reduced to two. The extensive algebraic calculations 
, that arise are done by a computer algebra system (CAS). The PDE's are then solved by a difference scheme using 
t— i ■ the Newton-Raphson algorithm of successive approximations for multidimensional equations. Additional features 
C*~) for global convergence are implemented. To obtain basic states that are sufficiently close to the solution, a one 
dimensional linearized version of the equation is solved with a numerical computation that reproduces the analytical 
results of Huang and Mura (1970). 

■ With this method, rigorous solutions of the nonlinear equations without any additional simplifications can be 
, obtained. The numerical results show a contraction of the singularity line which corresponds to the well-known 

^ Poynting effect in nonlinear elasticity. This combination of analytical and numerical computations proves to be a 
2? versatile method to solve nonlinear boundary value problems in complicated geometries. 



Ch ! 1 Introduction Bullough, and Smith (1955) that dislocation density is 

related to the differential geometric concept of torsion. 

H Topological defects in elastic media usually generate large After the result of Frank (1949) 1 the elastic energy of 

_ > . deformations which require finite elasticity for their de- topological defects in elastic continua has been consid- 

k> | scription. The theory of nonlinear elasticity goes back to e red in various publications on dislocations and discli- 

^ . Cauchy (1827), and significant contributions came from nations (Eshelby 1949; Kroupa 1960; Huang and Mura 

C3 Love (1927), Signorini (1930), and for the incompressible 1970; Nabarro 1979). Puntigam (1996) discussed topo- 

case, from Rivlin (1948a). For an introduction to nonlin- logical defects in a field theoretic context, 
ear elasticity, the reader may consider Beatty (1987) or However, calculations of the elastic energy of topolog- 

the standard textbooks Truesdell and Toupin (1960) and i ca i defects have always been restricted to the linearized 

Truesdell and Noll (1965). equations. In this case the problem can be treated using 

Topological defects, namely dislocations and discli- tensor analysis and stress function tensors of first and 

nations, have been investigated by means of continuum SCCO nd degree (Kroner 1959). Two assumptions usually 

theories (Kroner 1959; Kroner 1960) that were devel- — 7— — ; — ; , . , , 

' , 1 He showed that the energy increases relativistically with the 

oped after the discovery of Kondo (1952) and Bilby, dislocation velocity. 



1 



enter these approximations: 

Firstly, one restricts to linear elasticity in the sense 
that shearing stresses are assumed to cause a state of 
simple shear in the material (as in Huang and Mura 
1970). 

Secondly, the energy is assumed to be a function of 
the distorsion tensor that can be approximated by the 
gradient of the displacement field. Then, by applying 
Stoke's theorem, an integration over the boundaries only 
can be performed. Regarding the first point, however, it 
can be shown (e.g. Rivlin 1948c, p. 467), that in the 
general case shearing stresses alone cannot maintain a 
state of simple shear in the material. Also the second as- 
sumption is justified for small deformations only (Rivlin 
1948c, p. 476). In the region close to the defect core, 
where the deformations are large the linear approxima- 
tion fails to predict finite energies. Thus, these results 
are reasonable only outside the core region in where the 
energy density diverges. 

Large deformations in crystals have been investigated 
rarely (Frank 1951), and if, not by means of analytical 
but by general topological methods (Rogula 1976) with- 
out considering the elastic energy. 

There are several proposals how to treat nonlinear 
effects in elastic media with defects. Teodosiu (1982) ob- 
tained second-order effects for a straight screw disloca- 
tion by applying Willis' (1967; 1970) scheme which goes 
back to Signorini (1930). This requires, however, certain 
physical assumptions for every additional order of ap- 
proximation. Guo (2001) modelled numerically the large 
deformations of a hyperelastic membrane. The cylindri- 
cal symmetry allowed a the reduction to a onedimen- 
sional PDE. 

Lazar (2002a; 2002b) obtained solutions for edge and 
screw dislocations in an elastoplastic theory. 

Povstenko (1995; 2000) has treated the twist discli- 
nation with a nonlocal modulus. Even if his method in- 
volves numerical integration, it is based on an analytical 
derivation of the stresses, and does not allow a contrac- 
tion of the singularity line. 

Here a more general approach is proposed: The total 
elastic energy, i.e. the integral over the energy density 
must be a minimum under variation of the displacement 
field u. Additional constraints - like in the present case 
incompressibility - are included by means of a Lagrange 
multiplier. If the energy depends on first derivatives of 
u only, this leads to Euler-Lagrange equations of second 
order, even if the method given here allows the treatment 
of higher order equations. 

The general method outlined above is applied to a 
circular twist disclination (Huang and Mura 1970) in an 
incompressible hyperelastic continuum. It is a numeri- 
cally automatized process for obtaining rigorous solu- 
tions of the nonlinear equations without special physical 
assumptions. 

Section 2 gives a brief introduction to nonlinear con- 
tinuum mechanics, followed by a description of the circu- 



lar twist disclination. The appropriate coordinate system 
and the respective transformation of the the displace- 
ment vector is also explained there. All the numerical 
issues are addressed in section 3, and the results are dis- 
cussed in section 4. 

2 Analytical description 

2.1 Basic concepts of nonlinear contin- 
uum mechanics 

Nonlinear elasticity was founded in 1827 by Cauchy. 

One assumes an undeformed, euclidic continuum with 
Cartesian coordinates X = (X, Y, Z) (the so-called 'ref- 
erence configuration') and attaches in every point a dis- 
placement vector u that points to the coordinates x = 
(a;, y, z) of the deformed state ('configuration'): u = x 
X. ' 

From u one deduces a quantity that transforms the 
coordinates X of the undeformed state to those x of the 
deformed state: the deformation gradient 

F-= *i 
• 0X' 

or, in components, 

(1 + u x My u z \ 

v x 1 + vy v z , (1) 
w x v Y 1 + w z J 

where (u, v, w) denote the components of u and the sub- 
scripts differentiation. The symmetrical tensor 

B := FF T , (2) 

is called (right) Cauchy-Green tensor. It is convenient to 
introduce the so-called principal invariants (ii, I21I3) of 
a tensor that are defined as follows: 

h = Ai + A 2 + A 3 (trace) (3) 
h = AiA 2 + A 2 A 3 + A3A1 (4) 
h = AiA 2 A 3 (det) (5) 

(the Aj are the eigenvalues of B). Then, in general, the 
energy density W is a function of the principal invariants 

W = W{h,I 2 ,h) (6) 

of the Cauchy tensor B (e.g. Bcatty 1987), and the stress 
tensor T (defined as traction per surface element) is 
given by the constitutive equation 

T = /3 1 B + /3 E + /3_ 1 B- 1 (7) 

with the (3i{I\, I2, 13) being functions of the principal 
invariants. For small deformations, i.e. for Ji's close to 
1, and h-h close to 3, the /3's have a fixed values - that 
can be related to the known elastic constants in linear 
elasticity. 
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2.2 Incompressibility 

The nonlinear condition of incompressibility is given by 

h = A1A2A3 = det B = (det F) 2 = 1. (8) 

and not, as many texts on linear elasticity state, div u = 0. 
As a consequence, the elastic energy W depends on the 
principal invariants I\ and I 2 of B only and eqn. 6 re- 
duces to 

W = W{h,h) (9) 



The Euler-Lagrange equations are thus obtained from 



In the present paper, the special case 
W(h)=C(h-3) 



(10) 



is considered, which is called a Neo-Hookean material 2 
where C is a constant. 

2.3 Strain-Energy function expressed by 
displacements 

Considering again a Cartesian coordinate system and 
taking into account (1), (3) and (10), the energy per 
unit volume for an incompressible, nco-Hookean mate- 
rial reads 

W = C(e xx + e yy + e zz ), (11) 

where C corresponds to the shear modulus and the 
respective components of strain are given by (Cauchy 
1827; Love 1927, p. 60 and Rivlin 1948b, p. 461) 



e xx = u x + i(?i 2 +v 2 x + w 2 x ) 
£yy = v y + \{ u l + yl + wl) 

tzz = Wz + \{u\ + z\ + W 2 Z ). 



(12) 
(13) 
(14) 



Hereby u = (u, v,w) denotes the displacement vector in 
x-, y-, and ^-directions and u x , u y ... etc. the respective 
partial derivatives. 
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and 



d(x + u,y + v, z + w) 



1 = 0, 



(16) 
(17) 
(18) 

(19) 



d(x,y,z) 

which is equivalent to (8). Using the special form of W 
from (11) without a Lagrange multiplier still leads to a 
linear PDE. Adopting W from (11) is not an essential 
simplification of the problem, since the condition (19) 
is the main nonlincarity. The principles of the following 
calculations apply in the same way to a nonlinear W. 

2.5 Circular twist disclination 

A circular twist disclination is a nontrivial example of 
a topological defect which enforces large deformations 
of an elastic continuum. The twist disclination can be 
realized as follows (see fig. 1): 




2.4 Euler-Lagrange equations 

The energy of the topological defect can be found by 
minimizing W under the variation of the displacement 
field u which fulfill the correct spacial boundary condi- 
tions. In addition, the nonlinear constraint (8) is taken 
into account by means of a Lagrange multiplier A. This 
requires the minimization of the Lagrangian 

L = W + X(det F-l). (15) 

2 This corresponds to the condition ft = const. Note that in 
finite elasticity the energy density depends upon the functions /3;. 
Therefore, no general 'canonical' energy density exists. 

3 In an incompressible material, 3C equals Young's modulus E. 



Figure 1: Schematic description how to produce a twist 
disclination in an elastic continuum. The solid torus (size 
R, thickness r) is removed. Then the material is cut along 
the hatched surface. After twisting the cut faces by the 
amount of the Frank angle O, the material is rejoined by 
glueing. 

Imagine IF? filled with elastic material and remove 
a solid torus centered around the origin and with z as 
symmetry axis (fig. 1). Cut the material along the surface 
bounded by the inner circle with radius R—r of the torus 
in the x — y-plane (hatched disk in fig. 1). The cut faces 
are twisted with respect to each other by the Frank angle 
fl, and glued together again. 
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If one shrinks the removed torus to a singularity line, 
in the limit r — > one obtains a circular twist disclina- 
tion with radius R. In the subsequent computations, a 
sequence of finite values of r is used. 

The above description of the twist disclination is equiv- 
alent to the one given by Huang and Mura (1970), who 
used however the term 'edge disclination' 4 to indicate 
that the Frank vector is perpendicular to the disclination 
line, in analogy to edge dislocations where the Burgers 
vector is perpendicular to the defect line. They investi- 
gated analytically the twist disclination in linear approx- 
imation. 

2.6 Transformation to curvilinear coor- 
dinate systems 

The boundary conditions for a specific problem can be 
formulated most conveniently in an appropriately cho- 
sen curvilinear coordinate system. Furthermore, differ- 
ent scale factors in this coordinate systems allow to pave 
the regions of interest - for example regions of large de- 
formations - more densely with coordinate lines, which 
is of great numerical advantage. Frequently, one can use 
symmentries of a problem and reduce it to a lower di- 
mension. The variety of orthogonal coordinate systems 
allows to design an appropriate system for nearly ev- 
ery problem, even for those with less symmetry. To give 
some examples, straight line defects may bemodelled in 
a cylindrical system ( e.g. using a log-transformed ra- 
dial component), closed line defects with various elliptic 
coordinate systems. 

For the above described twist disclination, the toroidal 
system, as visualized in Fig. 2, is a suitable choice. To 
obtain the 3D-toroidal coordinate system from the 2-D 
bipolar system which is actually shown in Fig. 2, one 
must rotate the latter around the vertical axis. 

The transformation is given by 



cos^) sinh(?7) 

— cos(i9) + cosh(?y) 
sm{ip) sinh(ry) 

— cos(i9) + cosh(Ty) 

sin(tf) 

— cos(?9) + cosh(fy) 



(20) 
(21) 
(22) 



The parameter of the azimuthal rotation is 6 [0, 2n[, 
whereas the polar angle i9 ranges from to w. The hy- 
perbolic coordinate r\ ranges from — oo (left focus) to oo 
(right focus). 

4 To avoid confusion with screw and edge dislocations, the terms 
'twist' and 'wedge' are commonly used for disclinations (deWit 
1973a; 1973b; Zubov 1997). Rogula (1976) calls the circular twist 
disclination a 'third type of defect' and Unzicker (1996; 2000) a 
'screw dislocation loop'. In any case, the twist disclination re- 
ferred here causes locally a Volterra distortion of the 5th kind 
(sec Puntigam 1996). 



Toroidal Coordinate 
System 




Figure 2: Toroidal coordinate system obtained by ro- 
tating a bipolar coordinate system around the z-axis 
[r\ = 0). For symmetry reasons, only the region 77 > 
and < $ < 7r (at <j> — 0) is shown. The value i) = 00 
corresponds to a focus of the bipolar system or, after ro- 
tating around the z-axis, to a circular singularity line in 
the toroidal system. r\ = -d = is the infinitely far point. 

In the region in the vicinity of the circular singular- 
ity the Jacobian determinant is very small, rj = ■& = 
corresponds to the infinitely far point. The components 
of the displacement vector in 77-, 7?-, and ^-direction arc 
denoted as u, v and w. 

The transformation of differential operators like div 
and curl in curvilinear systems is long known (Love 1927, 
p. 54) and is easily done by CAS 5 . 

To transform the expressions for the energy W and 
for the condition (19) in curvilinear coordinates, how- 
ever, this is not sufficient, since the displacement vector 
components (u,v,w) are involved. One has to take into 
account that the basis vectors in a curvilinear system are 
subject to change and differentiate them. This technique 
is well-known in differential geometry since it is needed 
for a curved spacetime in the general theory of relativity. 

The spatial derivatives of the displacement vector 



u 



i l = (u,v,w) are . Taking into account the 
transformation of the basis vectors, the usual derivative 
has to be replaced by the 'covariant' derivative: 



D 
1 

Dx k 



d 

dx k 



(23) 



The connection T- l k is calculated from the Jacobian ma- 



trix A; (A, 1 = £ etc.) by 



dx m 3 



(24) 



5 Here the Mathematical^) package 'VectorAnalysis' has been 
used. 

6 In contrast to eqns. (16-18), (u, v, w) and refers to jy-, 1?-, 
and 0-directions in the toroidal system. 
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B/ is the inverse of A> (Schouten 1954, p. 121ff.).Eqn. (24 
has to be applied to all derivatives occurring in (16-19). 
This again is done by a CAS. 

The transformation considerably complicates the equa- 
tions. For example (11) becomes in toroidal coordinates 

W = u v (cosh n — cos i9) + v$ (cosh rj — cos #) + (25) 
u (csch rj — cos i? coth rj) — 2v sin t? — u sinh 77 + 
(tu 2 (cosh rj — cost?) 2 + i*4 (— cos t? + cosh??) 2 + 
w 2 (cos cosh 77 — l) 2 csch rf + w 2 sin -d 2 + 
(v v (cosh rj — cos i?) + u sin i?) 2 + 
(u v (cosh 77 — cost?) — 7; suit?) 2 + 
(w(csch 77 — (cos 1? coth 77)) — vsini3) 2 + 
(v$ (cosh 77 — cos 1?) — u sinh 77) 2 + 
(7j^(cosh77 — cost?) + usinliTi) 2 ) /2 

The transformation of cqns. (16-19), i.e. the full non- 
linear PDE to solve, is not written out here, 7 they amount 
to a text file of about 30 kB. The extensive algebraic cal- 
culations necessary to transform the equations are done 
by a CAS. 

2.7 Boundary conditions 

The circular singularity of the twist disclination corre- 
sponds to 77 = 00. A toroidal core 77 > rj max sourround- 
ing the singularity 77 = 00 is removed from the material. 
The chosen values of r] max = 2 . . . 3.625 correspond to 
core radii r = 0.23 . . . 0.049. Since 77 = i9 = is the in- 
finitely far point, only a region 77 > rj min and 7? > t? min 
is considered. 

The distance R of the singularity line to the symme- 
try axis 77 = (the size of the defect) is set to 1. 

The cylindrical symmetry of the problem in the co- 
ordinate ip reduces the 3D-problem to two dimensional 
one. Therefore, u, v and w depend on 77 and 7? only. Fur- 
thermore, due to mirror symmetry 7? is restricted to the 
region < 7? < ir. 

The boundary values of u, v and w are left free where 
ever possible in order to allow relaxing to the configura- 
tion of minimal energy. This is done at the surface of the 
torus sourrounding the singularity at 77 = T] max (torus 
shown in fig. 1), and at 77 = 0. The boundary conditions 
at t? — are enforced by the symmetry of the problem 
and at 7? = 7T by the 'cut-and-glue' condition with the 
Frank angle ft. 

Thus at the disk defined by 7? = n, v vanishes and u 
and w are fixed to a purely circular displacement corre- 
sponding to the Frank angle (see fig. 2). During en- 
ergy minimization, the shape of the removed torus is 
kept fixed, otherwise the material could overlap, which 
is physically impossible. The boundary values w($ = 

7 Note that minmization of W in the neo-Hookean case eqn. 11 
still leads to a linear PDE. Only the constraint of incomprcssibility 
(8) enforces the substantial nonlincarity of the final PDE. 



= Vmax) = u(i3 = tf min ,ri = i] max ) implicitly imple- 
ment this condition. Only for infinitesimal small values 
of the u component would vanish. In the linear approx- 
imation, 77, v and w vanish at 7? = 0. At the symmetry 
axis rj = the condition of incomprcssibility causes the 
vanishing of u and w without an explicit setting to zero. 

3 Numerical methods 

3.1 Discretization with a 
difference scheme 

The nonlinear Euler-Lagrange equations generated by 
the CAS are discretized on a twodimensional lattice by 
evaluating the coefficients rj, sin 7? etc. at every lattice 
point. The displacements and their derivatives are still 
maintained in analytic form. One obtains 4 variables 
(displacement vector u and Lagrange multiplier A) at 
each of the n grid points. Therefore 4 77 nonlinear equa- 
tions have to be solved simultaneously. 

3.2 Newton's method in multidimensions 

The solutions of the arising nonlinear system are ob- 
tained by the Newton-Raphson method of successive ap- 
proximations (Press 1993, chap. 9.7), which is a multidi- 
mensional extension of Newton's root finding algorithm. 
As in one dimension the algorithm starts from an ba- 
sic state and uses first derivatives to iteratively calculate 
approximations of the solution. 

The numerical computation of the An derivatives of 
the Euler-Lagrange equations is however numerically in- 
hibitive. This problem is solved by analytically differenti- 
ating the Euler-Lagrange equations at every lattice point 
with respect to the displacements u,v and 77; and their 
derivatives using a CAS. Thus the 18 quantities A^, A#, 

U v , U#, V v , V#, W. q , W#, U v $, V n s, Wrffi, U##, V m , W„, TTJtftf 

entering the Euler-Lagrange equations are discretized by 
difference molecules (Bronstcin 1985, p. 767f.). 

Only the simplest difference molecules are implemented 
to avoid instability of the Newton-Raphson method due 
to unrealistic smoothness assumptions. The possible grid 
dimensions do not allow to assume smoothness of the 
solution because the distorsions in the vicinity of the 
toroidal core region usually become very large and higher 
order difference operators contain no next-neighbor cou- 
pling. 

In addition, the difference scheme allows for a conve- 
nient formulation of the boundary conditions. 

3.3 Global convergence 

Although the above implementation is numerically very 
efficient, the global convergence of the direct multidimen- 
sional Newton-Raphson method is still not guaranteed. 
To obtain the global convergence of the algorithm one 
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usually defines a functional V on the solution space as 
the square of the l.h. sides of eqns. 16-19. If the New- 
ton step increases V with respect to the previous state, 
the algorithm used here 'walks back' along the Newton 
direction looking for a one-dimensional local minimum 
of V. (cfr. Press 1993, 9.7). The existence of such a lo- 
cal minimum is guaranteed, since at the basic state, V 
by definition of the gradient decreases along the Newton 
direction. It turns out that convergence improves if the 
'walk back' already is undertaken if V decreased slightly, 
and the Newton step was 'accepted' only if V decreases 
by less than a factor of 10. A factor 10 is easily reached in 
the vicinity of the solution, where the Newton-Raphson 
algorithm converges quadratically with distance. If the 
Newton step is accepted, its result is used as the basic 
state of the next iteration step around which the PDE is 
again linearized. In contrast to the advice given in Press 
(1993), it is found here that it is useful to minimize V 
along the Newton direcetion precisely. Although this re- 
quires more function evaluations, it performs superior to 
taking some premature value for a new linearization. 

For the one-dimensional minimization along the New- 
ton direction, the golden sectio algorithm (Press 1993, 
chap. 10.1) is used. Where the function is sufficiently 
flat 8 (functional V at three points vary by less than a 
factor 1.2), the golden sectio is accompanied by parabolic 
interpolation (Press 1993, chap. 10.2) of the minimum. 
This combination, is found to perform quite well in find- 
ing solutions. Moreover, it is not sensitive to small vari- 
ations of the above parameters. 

3.4 Basic states 

The success of the Newton-Raphson method critically 
depends on the basic state used for the first linearization. 
In the multidimensional case there is little hope that the 
algorithm directly converges to the solution without the 
above mentioned global-convergence methods. But even 
then the algorithm ends up in a local minimum if the 
first basic state is too far from the solution. There is 
for example, no chance to find a solution when starting 
with a nonspecialized basic state, e.g. with all variables 
set to zero. Sometimes linear approximations, e.g. states 
obtained by assuming small deformations, are chosen as 
basic states. In the present case, the linearized constraint 
of incompressibility div u = 0, still is too far from real- 
ity to be used for constructing a first basic state. Thus, 
the incomprcssibilty condition is directly implemented 
by using the symmetries of the toroidal coordinate sys- 
tem. In first approximation, it can be assumed that the 
displacement takes place along concentric circles with 
respect to the symmetry axis. Only for small displace- 
ments this corresponds to the tangential displacement 
of the linearized equation. Thus, no dilatation at all is 
allowed in the circular approximation. For practical rea- 

8 And if the point in the middle was not too excentric. 



sons, only the component w of the diplacement is var- 
ied. The components u and v are adjusted to meet the 
above requirement. Thisprocedure requires to solve only 
a onedimensional linear PDE for calculating the basic 
states. 

3.5 Interpolation and extrapolation 

For finer grids and smaller core radii, it becomes increas- 
ingly difficult to find solutions by means of the above 
method. 

It turns out that grid refinement by interpolating a 
previously found solution on a coarser grid is by far more 
efficient than starting with an independent new solution 
of the linearized system. The interpolation is performed 
by a 2D-spline algorithm (Press 1993, chap. 3.6) and the 
interpolated function evaluated at the refined grid - is 
taken as new basic state. Starting from this state, a new 
solution usually is obtained by Newton-Raphson. How- 
ever, this refined solution still has the same core radius 
as the previous coarse solution. For extending the grid 
towards smaller core radii, that is to larger rj maX} a very 
efficient method is to attach one lattice line to the old 
solution by linear extrapolation, and to use this extended 
solution again as a new basic state. All solutions for very 
small core radii have been obtained in this way. 

3.6 Matrix inversion and function eval- 
uation 

Matrix size increases very rapidly with the number of 
grid points. For example, the four functions (A, u, v, w) 
on a 30 x 30 lattice lead to a 3600 x 3600 matrix with 
almost 13 x 10 6 coefficients, which has to be inverted 
at every Newton step. Since the matrix is sparse, the 
inversion is more efficiently done by specialized sparse 
matrix algorithms. Here the corresponding packages of 
MatLab© are applied over a data exchange interface 
with Mathematica©. 

Since CAS are slow in evaluating trigonometric func- 
tions, the evaluation is done by an external C routine. 
This speeds up the computation of the matrix coefficients 
by a factor 10 as compared to the CAS. 

With this features, satisfactory results are obtained 
by running the program on a usual PC. 

4 Results 

4.1 Linear approximation 

In linear approximation the twist disclination can be 
treated analytically. Huang and Mura (1970) have found 
that the total elastic energy is 

W = tt 2 i? 3 ){[2 + (1 - Lf}K (1 - 1) - (26) 

2[1 + (1-1) 2 ]E(1-1)}, 
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Figure 3: Comparision of the numerical results for the 
linear approximation (filled circles) with the analytical 
results (solid line) of Huang and Mura (1970) in case of 
a core radius of 0.09 [r\ max = 3, r] min = tf min = 0.05, 
R = 1, \i — 3). The defect energy increases with the 
amount of the Frank angle f2 (ranging from to tt). For 
small fl, the results agree. 

where K and E denote the complete elliptic integral of 
the first and second kind. Eqn. 26 can be approximated 

W=^Vl 2 R 3 (±\og(8R/r)-^). (27) 

This result is used to test the algorithm for calculating 
the basic states (section 3.4). Fig. 3 shows that numerical 
and analytical results agree for small CI. The predicted 
fi 2 -dependence in eqn. 26 is reproduced by the linear 
algorithm. 

For larger values of fl, however, the energies calcu- 
lated numerically lie below the parabola of Huang and 
Mura (1970). This illustrates the two aspects of a linear 
approximation: both in the analytical and in the numer- 
ical treatment shear stresses are assumed to cause pure 
shear, but the latter approach lifts the small-deformation 
assumption by taking the energy function (11). 

4.2 Solutions of the complete nonlinear 
PDE 

The linear assumption of shear stresses causing shear 
deformations only is lifted now. The solutions of the fully 
nonlinear equation presented here are constraint to the 
the case Vl = tt where very large deformations occur. 
Fig. 4 visualizes these deformations for a core radius r = 
0.23. The other parameters are /i = 3, R = l,f] m i n — 
$min = 0.05. The dottet line indicates the singularity. 
Fig. 2 shows the corresponding undeformed state. 

The components of the displacement vector u, v and 
w, for the same parameters as in fig. 4, are shown in 
fig. 5. Note that the grid shown in fig. 5 is equidistant 
in the toroidal coordinates r] and i?, but not in a carte- 
sian frame, as it is evident from figs. 2 and 4. There, 




Figure 4: Solution of the nonlinear PDE. The respec- 
tive undistorted state is shown in fig. 2, with singu- 
larity as dottet line. A very large deformation with a 
Frank angle fl — tt is implemented by twisting the lower 
boundary $ = tt by the amount of tt/2. The lower part 
— tt/2 < & < (not shown here) is twisted in the oppo- 
site direction. The core radius is r about 0.23, f] max = 2, 

Vmin $min 0.05. 
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Figure 5: vector components u(rj, $), i>(?7,i?) and w(r), i?) 
for the solution fig. 4 (core radius r = 0.23). 
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Figure 6: Distribution of the elastic energy for an circular 
twist disclination with rj max = 2.5 (equivalent to r = 
0.145), and r) mm = 4m = 0.05. 

the lattice spacing is much finer in the region of large 
deformations. 

Even if it demands quite a lot of imagination, the 
reader should verify that the displacement components 
in fig. 5, attached to fig. 2, yield indeed the deformation 
fig. 4. 

In fig. 6 (r = 0.145) the distribution of the elastic 
energy plotted. At each grid point the energy density, 
weighted by the corresponding toroidal volume element, 
is shown. Most of the energy is stored in the region 
nearby the cut surface, where its density is greater than 
in the vicinity of the core region r\ — r\ max . This effect is 
even more pronounced for smaller core radiii r. 

4.3 Elastic energy 

Fig. 7 shows the total elastic energy of a circular twist 
disclination in function of r) max which parametrizes the 
core radius. 

Using the global minimization techniques discussed 
above, solutions over a wide range from r\ max = 2.0 
(which corresponds to r = 0.23) to i] max = 3.625 ( or 
r = 0.049) arc found. Fig. 7 shows the dependence of the 
energy (filled circles) from -q. 9 It is considerably reduced 
in comparison to the linear approximation of Huang and 
Mura (1970) (solid line). Although from numerical re- 
sults like those of fig. 7 it cannot be deduced that the 
energy approaches a finite value in the limit r — > 0, fig. 6 
suggests the energy density decreases sufficiently for con- 
vergence. 

The total elastic energies shown in fig. 7 are calcu- 
lated from a lattice of 16 points in ■&- direction, whereas 
the number of points in ry-direction ranged from 17 to 

9 r/ stands for r) ma x + Vmin here. 
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Figure 7: Dependence of the total elastic energy of cir- 
cular edge disclinations for rj max = 2.0 (r = 0.23) to 
Vmax — 3.625 (r = 0.049). The filled circles are obtained 
from the nonlinear modelling, the solid line shows the 
analytic result of Huang and Mura (1970). The cutoff of 
the infinitely far point was r\ m in = $min = 0.05, which 
corresponds to a distance of about 20 from the origin. 
There were 16 lattice points in -^-direction and a range 
from 17 to 30 for r\. 

30. The outer limit of the modelled region is defined by 
Vmin = fl-min = 0.05. Other solutions with a different 
number of grid points, and smaller values for r\ m i n and 
®min (0.01, 0.02) showed only a very slight dependence 
on these parameters. 



linear case 



v 




Figure 8: Example of a small core radius r = 0.09. Vec- 
tor component v in ^-direction for a solution of the lin- 
earized equations. 



4.4 Poynting effect 

The solutions of the nonlinear case show a characteristic 
phenomenon which only occurs in the elasticity theory 
of finite deformations, the so-called Poynting effect. 

The Poynting effect in its classic form is (Truesdell 
and Noll 1965, p. 193): 'When an incompressible cylin- 
der, free on its outer surface, is twisted, it experiences an 
elongation ultimately proportional to the square of the 
twist.' 

To illustrate the difference between the linear approx- 
imation and the nonlinear solution, fig. 8 and 9 show the 
vector component v in ^-direction for r = 0.09. In the 
region close to r\ = 0, i.e. in the vicinity of the symmetry 
axis, v has negative values and the displacement vector is 
pointing upwards in the region where z > 0. (cfr. fig. 2). 
By mirror symmetry it points downwards in the region 
of z < (— 7r<i9<0).Asa consequence the material is 
stretched and incomprcssibility forces the central region 
close to the cut surface to contract. 

The elongation indicated for z > is observed in the 
nonlinear case only. It is a displacement normal to the 
plane on which the stress is induced (the cut surface) . In 
the linear case from the assumption of pure shear caused 
v = at the symmetry axis (fig. 8). 

The Poynting contraction apparently approaches a 



nonlinear case 




Figure 9: As fig. 8, but for the nonlinear case. In the re- 
gion nearby i] — 0, i.e. in the vicinity of the symmetry 
axis, v has negative values, that means it is pointing 
upwards. The material is stretched and the cut surface 
contracts. 
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Figure 10: Poynting contraction as fraction of the orig- 
inal size of the defect. The contractions seems to reach 
a lower bound at .58, at a core radius of 0.049, which 
corresponds to rj max — 3.625. 

positive lower limit for very small core radii r. The results 
shown in fig. 10 suggest that the contraction does not fall 
below 58% of the original radius. 

4.5 Problems 

For both r) max < 2.0 and rj max > 3.75 it becomes increas- 
ingly difficult to find solutions of the nonlinear equations. 
In the first case the difficulty apcars to be related to 
the imposed boundary conditions which fix the torus di- 
ameter r (i.e. the core radius). In the second case the 
Poynting contraction leads to a net change in the torus 
volume, which apparently leads to numerical problems. 

On the other hand, for small core radii, the program 
finds solutions which satisfy the numerical convergence 
criterium, but show a physically unreasonable zigzag- 
pattern in the displacements. Possibly the occur in huge 
deformations (a volume element of the size of the core 
radius becomes stretched to the length of the singularity 
line) become numerically untractable at least with the 
present method ('relatively' simple coordinate system). 
However, there still may be a deeper physical reason for 
this instability. 

A simpler analogue to the circular edge disclination 
in case of large deformations is a twisted rubber band. 
Everyday experience shows that while the Poynting ef- 
fect only slightly elongates the band, the band sponta- 
neously looses its axial symmetry and coils up as soon 
as the twist surpasses a critical value. 

In case of the edge disclination, the present model - 
which essentially depends on rotational and mirror sym- 
metry is not able to correctly model such a transition. 
It is therefore not astonishing that meaningful solutions 
are found up to a maximal rj max only. 



5 Conclusions 

By combining analytical and numerical approaches, a 
general method for obtaining rigourous solutions of non- 
linear boundary value problems with nontrivial geome- 
tries has been developed. 

With this new method, for the first time the fully 
nonlinear problem of a circular twist disclination in an 
hyperelastic imcompressible material could be treated. 

While the extensive calculations necessary for formu- 
lating Eulcr-Lagrange-equations in a curvilinear coordi- 
nate system were done by the computer algebra system, 
the solution of the discretized equations was obtained by 
the Newton- Raphson method with additional global con- 
vergence features. The results give a nonlinear correction 
to the total clastic energy obtained in analytical treat- 
ments (Huang and Mura 1970). Moreover, the method 
allowed for the first time to model the contraction of 
the singularity line and to observe the Poynting effect, 
which is characteristic for for large deformations in finite 
elasticity. An application of the outlined method to to 
other nonlinear PDE's with boundary value problems in 
several dimensions seems possible. 

Acknowledgement. We are grateful to Dr. Markus 
Lazar for hints to the literature. 
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